Measuring Spatial Distribution of Local Elastic Modulus in Glasses 
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Glasses exhibit spatially inhomogeneous elastic properties, which can be investigated by measuring 
their elastic moduli at a local scale. Various methods to evaluate the local elastic modulus have 
been proposed. A first possibility is to measure the local stress-local strain curve and to obtain the 
local elastic modulus from the slope of the curve, or equivalently to use a local fluctuation formula. 
Another possible route is to assume an affine strain and to use the applied global strain instead of 
the local strain for the calculation of the local modulus. Most recently p], 0], Sollich and Barra 
introduced a third new approach which is easy to be implemented and has the advantage of low 
computational cost. In this contribution, we compare these three approaches by using the same 
model glass and reveal the differences among them caused by the non-affine component of the local 
strain. 

PACS numbers: 62.20.de, 62.25.-g, 71.55.Jv 
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I. INTRODUCTION 

It is well documented that glasses exhibit spatially in- 
homogeneous elastic properties, with coexistence of hard 
and soft domains when the elastic properties are mea- 
sured at a local (typically ten atomic sizes) scale 
It has been suggested that the local heterogeneity in the 
elastic properties is closely linked to several unusual prop- 
erties of glasses, which include low-temperature thermal 
properties [f|, an excess vibrational density of states, 
known as the "Boson peak" 0, H| , and anomalous acous- 
tic properties [9l-fllj]. Theoretical models [H, have 
also been proposed to relate the boson peak and the as- 
sociated thermal and acoustic anomalies to a randomly 
fluctuating shear modulus. Elastic heterogeneity is also 
reflected by the existence of a strong non affine char- 
acter in the elastic deformation of the material 
the displacement field at small scale is not obtained from 
the macroscopic strain, but the atoms undergo an ex- 
tra relaxation described as a non affine displacement, 
which has long range spatial correlations due to the elas- 
tic character of the problem [TH, . The scale £ of the 
elastic heterogeneities can be assessed by measuring the 
local elastic properties as a function of a coarse grain- 
ing size, and monitoring the convergence towards macro- 
scopic properties [5[. The elastic continuum approxima- 
tion for the acoustic excitations breaks down on a meso- 
scopic wavelength comparable to £, where a marked re- 
duction of the sound velocity and strong scattering were 
observed [2(J[2lj]. Moreover, the localized plastic events 
[22], [23[ that lead to glass plasticity are related to the elas- 
tic heterogeneities, with a strong correlation between the 
plastic rearrangements and the spatial map of the local 
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shear modulus [5]. The concept of elastic heterogene- 
ity appears therefore as a key to understand mechanical 
(elastic or plastic) properties of glasses, and the quantifi- 
cation of these heterogeneities is of primary interest in 
simulation studies of amorphous systems. 

Various methods have been proposed to evaluate the 
local elastic modulus of materials. Probably the most 
natural one is to measure the stress-strain curve at a lo- 
cal scale. The local elastic modulus is then calculated 
from the slope (first-derivative) of the local stress with 
respect to the local strain, in the same manner as the 
macroscopic modulus. As for the macroscopic elastic 
constants, this approach can be implemented either by 
using at a local scale the statistical mechanical formulae 
d, [24], Hl| that are obtained from linear response theory, 
or by applying an explicit deformation to derive the local 
stress-local strain relation directly [H,[26j]. Another, more 
approximate approach is to assume an "affine strain" : 
the applied global strain, instead of the local strain, is 
used for the calculation of the local elastic modulus. In 
this approach, the local modulus is calculated from the 
slope of the local stress-global strain curve, so that a sim- 
ple and quick stress measurement gives immediately the 
elastic constants. A previous study [3[ reported that the 
local modulus calculated by this approach is qualitatively 
consistent with that obtained from the local stress-local 
strain curve, although there are some differences at a 
quantitative level. Furthermore, most recently, Sollich 
and Barra proposed a new approach to evaluate local 
elastic properties [H, 0] ■ In this approach, the material is 
" frozen" except for the "target" local region. The frozen 
part undergoes an homogeneous affine deformation cor- 
responding to an imposed external strain, and is not al- 
lowed to relax non-affinely. The local elastic modulus is 
measured from the stress response in the target region to 
the applied global strain. This third approach is easy to 
implement and has the advantage of low computational 
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cost. In Q it was used to relate cavitation in a uniaxi- 
ally strained polymer to "weak spots" in the local bulk 
modulus. 

The elastic modulus tensor C is composed of three 
components: C = C B - C N + C K 0, Hit S3- Thc first 
term C B is the so-called Born term, which corresponds to 
the instantaneous elastic modulus under a uniform affinc 
deformation. The second C N is the component due to 
the non-affine deformation: The non-affine or internal 
motions of particles give rise to a decrease of the modu- 
lus. The third C K is the contribution from the kinetic 
energy to the modulus, which is much smaller than the 
other two terms and can be neglected for dense systems. 
For amorphous materials, the non-affine term C N is an 
important contribution, comparable in magnitude with 
the affine term C B [3, [H, Q3] ■ The three approaches to 
measure the local modulus described above evaluate this 
non-affine component C N in different ways. In previ- 
ous studies , all of the three different approaches ap- 
peared to provide reasonable and qualitatively consistent 
values of the local moduli. However, since these studies 
used different systems, the comparison remained at 
a qualitative level. In this contribution, we apply the 
three approaches on the same system and compare the 
corresponding results quantitatively. We point out the 
differences among them due to the non-affine component 
of the displacement field, which is discussed in details. 

The paper is arranged as follows. In Sec. [Ill we briefly 
summarize the basic definitions of local bulk and shear 
moduli. In Sec. IIII1 our Lennard- Jones (LJ) glass system 
is described and the three approaches to measure the lo- 
cal modulus -fully local, affine strain, and frozen matrix- 
are presented. In Sec. IIV| we discuss and compare the 
results of the local moduli obtained by using the three 
approaches. Finally, we conclude with a summary of our 
findings in Sec. [V] 



W for our glass system (a 3-dimensional LJ glass). Ref- 
erence H also demonstrated for a 2-dimensional LJ glass 
that the symmetry of the local modulus tensor breaks 
down as the linear size W becomes small. 

In a three dimensional system there are six mutually- 
independent deformations represented by the following 
strain tensors 1301: 
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where e™, e™, e™2, e^>, e^|, and are the strengths of 
the strains. The strain tensor e™ represents an isotropic 
bulk compression of the domain m. Thc other five tensors 
express volume conserving shear deformations. The two 
tensors, e^j; and ej|, are the "pure" shear deformations 
(plane strain and triaxial). The three tensors, e^, e^, 
and e^, are the "simple" shear deformations. All the 
linear deformations are realized as superpositions of these 
six deformations (one bulk and five shear deformations) 

M- 

From the six deformations expressed in Eq. ([2]) , one de- 
fines six moduli: one bulk modulus and five shear moduli 
[30| . The bulk modulus K m is defined from the pressure- 
volume change relation under the isotropic bulk deforma- 
tion el n : 



II. LOCAL BULK AND SHEAR MODULI 

The local elastic moduli are defined in a manner similar 
to the one used for macroscopic moduli. Let us focus on 
a local cubic domain m of linear size W in a glass sample. 
The local elastic modulus in the cube m can be defined 
as the first-derivative of the local stress with respect to 
the local strain. The stress-strain relation in the domain 
m is written as 
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Here the pressure p m is the trace of the stress tensor, 
(cr™ + cr™ + cr™)/3, and the volume change SV m 
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3e™ In addi- 



is written as SV m /V m = e% 

tion, five shear moduli, G™, Gf,Gf, Gf, and Gg% are 
defined from the shear stress-shear strain relations un- 
der two pure shear deformations, and ej|, and three 
(1) simple shear deformations, e£|, e^, and e^, respectively: 



where <x m , e m , and C m are respectively the local stress, 
local strain, and local modulus tensors. The initial stress 
<r 0m generally has some non-zero values (28|. All the 
quantities in Eq. (p} are local ones, which depend on the 
position r and the size W of the cube m. Additionally, 
we remark that the local modulus tensor CKy is not nec- 
essarily symmetric with respect to the exchange of ij and 
kl like the macroscopic modulus Cijki [2{|. Indeed, we 
observed that G™ fc/ is not symmetric in cubes with small 
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where the five shear stresses are written as a™ = (cr™ 
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cr™ , respectively. These bulk and shear moduli 
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can be expressed as linear combinations of C, 
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the matrix C 



of 



to the diagonal components 
= p<"- l c m P m with P m = 

[ e T' e ^ e ^ e ^3> e % e ^\- In the case of isotropic sys- 
tems, like a LJ glass, the macroscopic modulus tensor 
C is characterized by one bulk modulus K = \ + 2p/'i 
and one shear modulus G = p 311, where A and p are 
the Lame constants. As the size W of the cube m be- 
comes large, the local bulk modulus K rn tends to the 
macroscopic value of K, and the five shear moduli, G™, 
G™, G™, G™, and G™, which generally have different 
values from each other, converge to the same value of the 
macroscopic G. 



III. METHODS 

A. System preparation 

The system considered in the present study is a 3- 
dimensional LJ monatomic glass model, described in Ref. 
[2lj . The interaction energy between two particles is 
(f)(r) = 4e[(er/r) 12 — (c/r) 6 ], where r is the distance be- 
tween the two particles, e is the depth of the potential 
well, and a is the particle diameter. The potential was 
cut-off and shifted to zero at r = r c = 2.5a. The sys- 
tem contains N = 4, 000 particles in a simulation box of 
constant volume V under periodic boundary conditions. 
In the following, all numerical values are expressed in 
LJ units: length in a, temperature in e/ks (&b is the 
Boltzmann constant), and time in r = (ma 2 / e) 1 / 2 (m is 
the mass of a particle). The number density was fixed at 
p = N/V = 1.015, which implies a linear dimension of the 
simulation box L — V 1 ^ = 15.8. For reference, we note 
that at the number density considered here, p = 1.015, 
the melting temperature of the system is T m ~ 1, and 
the glass transition temperature is T g ~ 0.4 32]. The 
glass phase of the system was realized as follows. We 
first equilibrated the system at the temperature T = 2 in 
the normal liquid phase by using a standard NVT molec- 
ular dynamics (MD) simulation. Next, we quenched the 
system down to T ~ 10~ 3 in the glass phase with a fast 
quench rate dT/dt = 4 x 10 2 . After quenching, the sys- 
tem was relaxed for a sufficiently long time to stabilize 
the total energy. 



B. Three approaches to measure local elastic 
modulus 



We subdivided the cubic glass sample (box size L = 
15.8) into a grid of size M x M x M and measured the 
local elastic modulus for small cubic domains (size W) 
centered at these grid points. In the present study, we 
set the cube size to be W = 3.16(= L/5), 5.27(= L/3), 
and 7.90(= L/2), and the grid size is set to be 40 x 40 x 
40 for W = 3.16 and 20 x 20 x 20 for W = 5.27 and 
7.90. The small cubes are distinguished by the index m 
(m = 1, 2, M 3 ). In the following, we describe the three 
approaches used to measure the local modulus. 



1. Fully local approach 

In this approach, one considers the local stress cr m and 
the local strain e m . The local modulus C m is obtained 
from the first-derivative of er™ with respect to e m , as in 
Eq. ([1]) . The approach can be implemented either by us- 
ing the equilibrium fluctuation formula [3j, [2J, |25j or by 
performing an explicit deformation to obtain the local 
stress-local strain relation directly 0, [2(| • In the present 
study, the fluctuation formula was used, as described be- 
low. 

Fluctuation formulae are obtained within the frame- 
work of equilibrium statistical mechanics [3, H3, HH, I27I 
l33l - [3"5l | . The local stress er™ for the small cube m is cal- 
culated as 
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where p m is the number density in the cube m, rf de- 
notes the vector joining particles a and b, and r ab is the 
distance between these two particles. The quantity q ab 
represents the length of the line segment rf 3 located in- 



If the vector 



is located outside the 

b l„ab 



side the cube m 

cube m, then q ab — 0. The term q ab / r ab determines 
the contribution of each pairwise interaction to the local 
stress crj™. The summation of crj™ over the entire system 
yields the usual macroscopic stress tensor a^ [36J: 
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The local modulus tensor C*" kl for the cube m is cal 



(a) Affine strain approach 



(b) Frozen matrix approach 



FIG. 1. (Color online) Schematic illustration of the simple shear deformation: (a) affine strain approach and (b) frozen matrix 
approach. The cubic box drawn by the yellow lines indicates the local cube m. In the affine strain approach (a), all the particles 
(red particles) are allowed to move non affinely. In the frozen matrix approach (b), only the particles in the local cube m (red 
particles) can move freely, whereas the particles in the other frozen part (blue particles) are restricted to only move affinely. 



culated from the following equations: 



2. Affine strain approach 
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the non-affine component 



contribution Cf-™i 
Cijki> calculated through the correlation between a^- 



(8) 

where () represents the ensemble average. As described 
in the introduction, there are three contribution to the 
elastic constants: the affine component (the Born term) 

and the kinetic 
We note that the non-affine term 

and 

<7ki, is not symmetric with respect to the exchange of ij 
and kl, and so the total modulus Cj™ kl is not a symmet- 
ric tensor in general. Earlier work |3j measured the local 
shear modulus of a polymer glass by using this method 
and revealed the inhomogeneous distribution of the mod- 
ulus. 

In order to apply Eqs. ©-([SI), the averages were per- 
formed on trajectories generated by standard NVT MD 
simulation on the glass sample at the low temperature 
T = 10~ 3 . The time step of the MD simulation was 
St = 5 x 1CP 3 , and the total duration of the runs was 
t = 10 5 (2 x 10 7 steps). The averages were performed 
over 10 4 configurations, separated by a time lag of 10 LJ 
units. 



In this second approach, one assumes an "affine 
strain", i.e., that the entire glass sample is strained uni- 
formly. This assumption says that the local strains e m 
of all the small cubes m are represented by the global 
strain e applied to the sample. The local modulus C m 
is defined and calculated based on Eqs. (H])-©, with the 
local e m replaced by the global e. This approach, which 
is obviously very simple to implement, ignores the spatial 
variations of the local strain field or, equivalently, its non 
affine component. 

To implement this approach we deform the entire glass 
sample in six ways: one bulk deformation and five shear 
deformations (two pure shear and three simple shear de- 
formations), which correspond to the six strain tensors 
in Eq. © with e m replaced by e. These deformations 
are illustrated schematically in Fig. [T] (a). For each de- 
formation, we calculate the local stress <r"f of the cube 
to as a function of the applied global strain e. Here, we 
use a formulation for <rj™ slightly different from that of 
Eq. ^ used in the fully local approach: a™j is calculated 
from the summation of the atomic stresses over the cubic 
domain m [3?], HH , 
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where the summation of a is performed over particles in 
the cube to. The same formulation of cr™ was used in 
Refs. 0,111. We note that the definition Eq. © of cr" 1 
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FIG. 2. (a) Pressure and (b) potential energy per particle versus strain 3et,, under the quasi-static bulk deformation €b- We show 
two deformation curves: afEne deformation with relaxation (solid curve) and affine deformation without relaxation (dashed 
curve). The macroscopic bulk modulus obtained from the slope of the stress-strain curve is K = 59.7 (AfEne + Relaxation) 
and K B = 60.2 (Affine). We also obtained the same values of K and K B from a quadratic fit of the potential energy. 



corresponds to the Eq. © with the term q ab /r ab 

1 (a,b G m), 

1/2 (a G 77i, b </ m), 

1/2 (a f m, b G m), 

(a, b £ m). 
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The usual macroscopic stress tensor atj is still recovered 
by the summation of a^- over the entire system, as in Eq. 
([7]) . The bulk modulus K rn is now obtained from the lo- 
cal pressure-global volume change relation, i.e., Eq. ([3]) 
with the local e™ replaced by the global t b . The shear 
moduli, G\ n , Gf, Gf, Gf, and Gg 1 , are calculated from 
the local shear stress-global shear strain relations, i.e., 
Eq. dU) with the local e™ replaced by the global e s . Ac- 
cording to the study [|[ j the assumption of "affine strain" 
can be qualitatively acceptable but causes quantitative 
deviations from the fully local approach. We examine in 
details the validity of this second approach in the follow- 
ing. 

In the present study, we have deformed the glass sam- 
ple by using a "quasi-static" protocol, i.e., the system 
was deformed at zero temperature T = 0, as described 
in Refs. [!, [H, H3- To achieve this, the system was first 
quenched using a steepest descent method from T = 1CP 3 
to T = 0, into the nearest energy minimum. Next, the 
system was submitted to an imposed deformation by ap- 
plying strain steps Se = 10 -5 with Lees-Edwards peri- 
odic boundary conditions [39j]. After each strain step 5e 
was imposed, the entire system was relaxed into its new 
closest energy minimum by the steepest descent method. 
The quasi-static deformation was performed until the ap- 
plied strain reached e = 0.002 (i.e., 0.2%). We note that 
e = 0.002 is small enough that the system deforms elas- 
tically, and the stress-strain curve is linear for all of the 
studied six deformations. We also remark that the differ- 
ence of the temperature between T — of this approach 



and T — 10 -3 of the fully local approach is very small so 
that this temperature difference is not expected to cause 
a discrepancy between those two calculations. 



3. Frozen matrix approach 

The third approach was originally introduced by Sol- 
lich and Barra [3,0. This also is an explicit type method, 
similar to the affine strain approach. We first "freeze" 
the system except for the "target" local region, i.e., the 
local cube m. The frozen region is restricted to deform 
only affinely due to the external strains, not allowed to 
relax non-affinely, whereas the target region m deforms 
non- affinely. Then we apply six types of deformations, 
one bulk and five shear deformations, on the entire sam- 
ple. We show a schematic illustration of this approach in 
Fig.[IJb), where the red particles are particles in the cube 
m, and the blue particles are particles in the frozen re- 
gion. During the deformations, the particles in the cube 
m can move freely, whereas the particles in the frozen 
region are displaced only affinely by the external strains. 
In this situation, where the parts surrounding the local 
cube m are frozen, the local strain e m of the cube m 
coincides exactly with the applied global strain e. For 
each deformation, we calculate the local stress a 7 ^ of the 
cube m as the function of the local strain e m = e. We 
note that the formulation of Eq. ^ (i.e., the summa- 
tion of the atomic stresses) is used for the calculation of 
crj™. Then, one bulk modulus K m and five shear moduli, 
G™, G™, G™, G^, and G£\ are determined from Eqs. 
([3]) and (|4]) . This approach is also easily mplemented, as 
the local strain e m = e is an input quantity. A recent 
numerical work [3] investigated the local bulk modulus 
of a glassy polymer by using this approach and obtained 
reasonable results in relation to cavitation events. 

We used a "quasi-static" protocol, as in the affine 
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FIG. 3. (a) Shear stress and (b) potential energy per particle versus strain 2e s under the quasi-static shear deformation e S 3. 
We show two deformation curves: affine deformation with relaxation (solid curve) and affine deformation without relaxation 
(dashed curve). The macroscopic shear modulus obtained from the slope of the stress-strain curve is G = 14.9 (Affine + 
Relaxation) and G B = 36.5 (Affine). We also obtained the same values of G and G B from the quadratic fitting of the potential 
energy. Note that we obtained five values of G ~ 15 and G B ~ 36 from the five shear deformations e s i, e S 2, e S 3, e S 4, and e s s. 
The values of G = 14.9 and G B = 36.5 are the averages over these values. 



strain approach. We first quenched the glass sample from 
T = 10~ 3 to T = using a steepest descent method, and 
then the quenched system was frozen except for the local 
cube to. The frozen system was submitted to an imposed 
deformation by applying strain steps 5e = 10 -5 , where all 
the particles move affinely. After each strain step Se was 
imposed, the system was relaxed by the steepest descent 
method. During the relaxation (energy minimization), 
only the particles in the target cube to (red particles in 
Fig. [1] (b)) can move toward the minimum energy point, 
whereas the particles in the frozen part (blue particles in 
Fig. Q] (b)) are stuck. The deformation was performed 
until the applied strain reaches e = e m = 0.002 (i.e., 
0.2%). 

IV. RESULTS 

A. Macroscopic elastic modulus 

We first investigated the macroscopic stress-strain re- 
lation and the macroscopic modulus, which are obtained 
from the quasi-static deformation. In Figs. [2] and [31 we 
plot the macroscopic stress as the function of the applied 
global strain: Fig. [2]for the pressure p under the isotropic 
bulk deformation e^, and Fig. [3] for the shear stress a S 3 
under the simple shear deformation e S 3 . The values of the 
bulk modulus K and the shear modulus G, which were 
obtained from the slopes of the curves, are K = 59.7 and 
G = 14.9. Note that we obtained five values of G ~ 15 
from the five shear deformations e s i, e S 2, e S 3, Cs4, and 
e S 5. The value G = 14.9 is the average over these five val- 
ues. The values of K = 59.7 and G — 14.9 are consistent 
with the previous work plj . where the longitudinal and 
transverse sound speeds, cl = \]{K + 4G/3)//0 ~ 8.8 



and ct = \J G/p ~ 3.8 (p = 1.015 is the mass density), 
were obtained for the same glass model as ours. In addi- 
tion, in Figs. [2] and we also plot the potential energy 
per particle. The potential energy is a quadratic function 
of the strain, and the moduli of K and G can be defined 
from the second-derivative of the potential energy p8l |. 
We obtained the same values of K ~ 60 and G ~ 15 
from the quadratic fitting of the potential energy. 



Also, Figs. [2] and [3] display the stress and the po- 
tential energy for the affine deformation, where the re- 
laxation (energy minimization) is not performed. For 
the isotropic bulk deformation in Fig. [21 the pressure 
profile is almost the same as that with relaxation. The 
bulk modulus K B = 60.2 (Born term) of the affine de- 
formation is nearly equal to K — 59.7. The potential 
energy also shows the same profile in both cases, with 
and without relaxation, respectively. This result indi- 
cates that the non-affine component is very small, and 
the bulk modulus is dominantly determined by the affine 
component. A previous study [13 obtained similar re- 
sults for a slightly polydisperse LJ glass. On the other 
hand, it is clearly observed that the relaxation causes 
a marked decrease of the shear stress and the potential 
energy under the shear deformation. The shear modu- 
lus G B = 36.5 (Born term) of the affine deformation is 
much higher than G = 14.9. The non-affine component 
G N = G B — G = 21.6 is of the same order magnitude as 
the affine component [3, EH, G3] ■ Therefore, the shear 
modulus has a large non-affine component, which is im- 
portant for amorphous materials. 
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FIG. 4. Distributions of bulk modulus and shear modulus for different cube sizes W = 3.16, 5.27, and 7.90 calculated by 
the three approaches described in the text: fully local approach ((a),(b)), affine strain approach ((c), (d)), and frozen matrix 
approach ((e), (f)). The shear modulus distribution P(G m ) is obtained by averaging the five distributions P(G?i I ), P(G™), 
P(G™), P(G™), and P(G™). We also show the Gaussian distribution functions fitted to each distribution (solid lines). The 
vertical solid lines indicate the values of the macroscopic moduli obtained from Figs. [2] and El i.e., K = 59.7 and G = 14.9. 
In (a)-(e), the average value is independent of the cube size W and coincides with the macroscopic value, whereas in (f), the 
average value varies with W and seems to tend to the macroscopic value with increasing W . 
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FIG. 5. Comparisons of the probability distributions for (a) bulk modulus and (b) shear modulus for the three approaches: 
fully local approach (white circle), affine strain approach (square), and frozen matrix approach (triangle). The local cube size 
is W = 3.16. The distribution of the Born term, calculated by the fluctuation formula Eq. (J5}, is also plotted for comparison 
(black circle). The vertical lines indicate the macroscopic values, K — 59.7 (solid line) and K B = 60.2 (dashed line) in (a), and 
G = 14.9 (solid line) and G B = 36.5 (dashed line) in (b). 



B. Distribution of local elastic moduli 

We next investigated the local moduli quantified by 
the three approaches described in Sec. IIII1 In Fig. [H we 
show the distributions of the local bulk modulus K m and 
the local shear modulus G m . The considered cube sizes 
are W = 3.16, 5.27, and 7.90. We confirmed that the five 
shear moduli, Gf , GV? , Gf , Gf, and Gf exhibit almost 
identical distributions, therefore in Fig. 2] we plot data 
averaged over these five components. The distributions 
calculated by the three approaches are all well fitted by 
gaussian distributions 04a] (solid lines). In Fig. |H we 
indicate the values of the macroscopic moduli obtained 
from Figs. [2] and 1 (i.e., K ~ 59.7 and G ~ 14.9) by the 
vertical lines. Except for the shear modulus distribution 
of the frozen matrix approach, all distributions exhibit 
an average value independent of the cube size W, and 
this average value coincides with the macroscopic one. 
The shear modulus distributions obtained in the frozen 
matrix approach, in contrast, exhibits an average value 
that depends on the cube size W, and seems to converge 
to the macroscopic value as W increases, although the 
convergence is rather slow. 

In Fig. [5l we show the comparison of the distribu- 
tions obtained from the three approaches for the cube 
size W — 3.16. In the same figure, we also plot the dis- 
tribution of the Born term (the affine component). Here 
it has to be noted that the Born term can be obtained 
either by the fluctuation formula through Cf^ in Eq. 
(|8]l or by the explicit way, i.e., performing explicit affine 
deformations as we do in Figs. [2] and [3] ("without re- 
laxation" case). We confirmed that these two methods 
produce identical distributions of the Born term. The 
distribution obtained by the fluctuation formula is also 
shown. From Fig. [SJa), it is evident that the three ap- 



proaches produce nearly identical bulk modulus distribu- 
tions, and these distributions coincide well with the Born 
term distribution. This result indicates that the non- 
affine component of the local bulk modulus is very small, 
which leads to the coincidence of the three approaches. 
As shown in Fig. [2] the macroscopic bulk modulus has 
a very small non-affine component and is mostly deter- 
mined by the Born term. Such situation is also true for 
the local bulk modulus. 

In contrast, the situation is totally different in the case 
of the local shear modulus. From Fig. [SJb), we clearly 
observe that the shear modulus distribution exhibits 
large differences among the three approaches. Also, the 
distributions are very different compared to the Born con- 
tribution alone, which means that the local shear modu- 
lus has a large non-affine component, as does the macro- 
scopic shear modulus shown in Fig. [3] The differences 
among the three approaches are indeed caused by this 
large non-affine component. In addition, there are two re- 
markable differences observed in Fig. [Sfb). Firstly, when 
comparing the distributions calculated by the fully local 
and the affine strain approaches, we see that the latter 
exhibits a narrower distribution (smaller variance). This 
result indicates that the spatial variations of the local 
strain field make the shear modulus distribution wider 
(more heterogeneous). Secondly, the frozen matrix ap- 
proach exhibits a much larger average value than the 
other two approaches, whose average values both coin- 
cide with the macroscopic one. We believe that the larger 
average value is caused by the additional constraint re- 
sulting from freezing the environment, which limits the 
non-affine motions of particles in the cube m. As the 
cube to can not be fully relaxed during the energy mini- 
mization, a larger stress and shear modulus are obtained. 
To support our explanation, we also observe that the av- 
erage value of the frozen matrix approach lies between 
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FIG. 6. Comparison of average and variance of the distributions calculated by the three approaches and of the Born term: 
fully local approach (open circle), afnne strain approach (square), frozen matrix approach (triangle), and the Born term (closed 
circle). The averages are shown in (a) and (b), and the variances are presented in (c) and (d). The horizontal lines indicate 
the macroscopic values, K — 59.7 (solid line), and K B — 60.2 (dashed line) in (a), and G = 14.9 (solid line) and G B = 36.5 
(dashed line) in (b). 



the macroscopic values G = 14.9 and G — 36.5, which 
correspond to the values of the non-constrained system 
and the fully-constrained system, respectively. This re- 
sult is clearly consistent with the interpretation that the 
local cube m is only partially relaxed. 

In Fig. |51 we compare averages and variances of the 
distributions calculated by the three approaches for three 
cube sizes, W — 3.16, 5.27, and 7.90, respectively. From 
this figure, we can emphasize the differences among the 
three methods more quantitatively. For the local bulk 
modulus, the three approaches exhibit nearly same aver- 
age and variance values for all three Ws, and these values 
coincide well with those estimated from the Born term 
only. The average values agree with the macroscopic val- 
ues K and K B (K ~ K B ). On the other hand, in the 
case of the local shear modulus, the affine strain approach 
shows a smaller variance value compared to the fully local 
one. Again, this is because affine strain approach does 
not consider the spatially varying local strain field. In ad- 
dition, the frozen matrix approach exhibits much higher 
average values than the other two approaches. The aver- 
age value of the frozen matrix approach lies between G 
and G B and seems to converge to the value G asW gets 
large. The constraints induced by the freezing of the ma- 
trix prevents the non affine field from fully contributing 



to the elastic constants, and causes such large average 
values. 

C. Spatial distribution of local elastic modulus 

We have also compared the spatial distributions of the 
local modulus, represented by color maps, for the three 
approaches. The spatial maps are shown for the case 
W = 3.16 in Fig. [7J Here, we visualize the maps for the 
variables K m and G™ , which are normalized by both the 
averages and the variances: 



var (11) 
G™ = r , (i€ 1,2,3,4,5), 

where the subscripts "ave" and "var" mean average and 
variance, respectively. By considering normalized vari- 
ables K m and G™, we are able to emphasize domains of 
the system which are relatively soft or hard. In the same 
figure we also show the spatial map of the Born term 
alone, which is also normalized as in Eq. (|11[) . 

For the local bulk modulus K m , the spatial maps of 
the three approaches correlate well with each other and 
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FIG. 7. (Color online) xy-layers taken from the spatial distributions of the bulk modulus K m and the shear moduli G^ 1 , Gj 1 
obtained by three approaches: fully local approach ((a)-(c)), affine strain approach ((d)-(f)), and frozen matrix approach ((g)- 
(i)). The value of the local modulus is normalized by its average and variance, i.e., X m = (X m — X avo )/X var (X £ K,Gi, G3) 
(see Eq. We also show the spatial maps of the Born terms ((j)-(l)), for comparison. X and Y coordinates are presented 

in units of the box length L — 15.8a. The Z coordinates of all the layers correspond to L/2. 
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FIG. 8. Comparison of the distributions of (a) bulk modulus and (b) shear modulus for two systems with N = 4, 000 (open 
circles) and N = 32,000 (closed circles) particles. The distributions were calculated by the fluctuation formula Eq. © (i.e., 
fully local approach). The local cube size is W = 3.16. 



are all very similar to that of the Born term K Bm . Again, 
this result indicates that there is a very small non-affine 
component in the local bulk modulus, which leads to the 
very small differences among the three approaches. On 
the other hand, differences among the three approaches 
are observed in the maps for the shear moduli G™ (pure 
shear) and G™ (simple shear). The large non-affine com- 
ponents of the shear moduli indeed change the soft/hard 
character of different domains in the spatial maps of the 
Born terms Gf m and Gf m . The fully local approach 
and the affine strain exhibit rather similar soft and hard 
parts, whereas the spatial map of frozen matrix approach 
looks like that of the Born term, rather than those of the 
other two approaches. This observation also indicates 
that the non-affine component of the frozen matrix ap- 
proach is limited, due to the constraint induced by the 
affinely deformed matrix, so that the modifications with 
respect to the Born term map are limited. 

From the above results (comparison of probability dis- 
tributions and spatial maps), we conclude that all the 
three approaches can be safely used to measure the local 
bulk modulus, which is characterized by a small non- 
affine component. In contrast, in the case of the local 
shear modulus, which has a large non-affine component, 
it is more appropriate to use a fully local approach to 
deal with both the local stress and the local strain and 
measure the non-affine component correctly, without any 
constraint. 

We conclude this section with two remarks. First, our 
results are qualitatively consistent withprevious studies 
0, 0] on different systems. The study [S[ observed that 
the values of the local shear modulus obtained from the 
local stress-global strain relation (i.e., affine strain ap- 
proach) are qualitatively consistent with those obtained 
from the fluctuation formula Eq. © (i.e., fully local). In 
our study, the fully local and the affine strain approaches, 
show similar modulus distributions. Indeed, although the 



two approaches exhibit modulus distributions which are 
different at a quantitative level, the differences are not 
so large, as evident from Fig. [SJ Moreover, the spatial 
maps corresponding to the two methods are very simi- 
lar (see Fig. [7]). This observation is analogous to what 
noticed in [3j]. In addition, the study [J] obtained very 
convincing results by using the frozen matrix approach 
but only focused on the local bulk modulus, ignoring lo- 
cal shear modulus. Therefore, our present results are 
consistent with and substantially expand previous simi- 
lar works [1, H| • 

Second, we briefly discuss now the impact of system 
size on the calculation of local moduli. Besides the sys- 
tem formed by N = 4, 000 particles which we have used 
in this study, we also considered a larger system, with 
N = 32,000, to check for system size effects. We calcu- 
lated the local moduli for this larger system by using the 
fluctuation formula Eq. (JSJ (fully local approach). We 
compare the modulus distributions for the two sizes in 
Fig. [5] It is evident that the two systems exhibit the 
same distributions for both the bulk modulus and the 
shear modulus. This result indicates that there are no 
system size effects on the local modulus. However, the 
total length of the MD trajectory used for performing 
the ensemble average of Eq. © has been found to be 
very different in the two cases. While a length t = 10 5 
was sufficient to obtain converged results for the small 
system of N = 4, 000, a length of t — 10 6 was required 
for the larger system of N = 32,000. More specifically, 
the correlation term between local and global stress, 
(cr™CT/c;), in the non-affine component G^j 1 , necessitates 
a longer sampling time to be estimated correctly. Indeed, 
the larger system is characterized by short wavelength 
modes, which contribute to fluctuations of the local stress 
a?? and evolve on slow time scales, which must be cor- 
rectly sampled. Concluding, convergence with simulation 
time must be carefully checked for larger systems. 
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V. CONCLUSIONS 

In the present study, we have applied three different ap- 
proaches, "fully local" , "affine strain " , and "frozen ma- 
trix" , to measure the local elastic moduli, bulk modulus 
and shear modulus of a LJ monatomic model glass. For 
the case of the local bulk modulus, the three approaches 
give nearly identical probability distributions and spatial 
maps. This is because the non-affine component in the 
bulk modulus is relatively small. The value of the bulk 
modulus is mostly determined by the Born term (the 
affine component), and therefore only small differences 
among the three approaches are observed. However, the 
situation becomes notably different for the local shear 
modulus. The three approaches exhibit different distri- 
butions and different spatial maps, even at a qualitative 
level. In the case of the shear modulus, the non-affine 
component has the same order of magnitude as the affine 
component, which causes the large differences among the 
three approaches. 

The difference between the fully local and the affine 
strain methods comes from the use of the local strain or 
global strain for the calculation of the local modulus. In 
the affine strain approach, the use of the global strain in- 
stead of the local value results in a variance narrower than 
that obtained from fully local approach, where the spatial 
variations of the local strain field is taken into account. 
In addition, the difference between the frozen matrix ap- 
proach and the other two approaches arises from that 
the system, except for the target local cube, is frozen. In 
the frozen matrix approach, the constraint of the affine 
displacement applied to the matrix severely restricts the 



relaxation in the local region, and the non affine compo- 
nent of the shear modulus is therefore underestimated. 
As a result, the average value of the shear modulus is sig- 
nificantly larger than in the other two approaches. The 
spatial map of the shear modulus obtained in the frozen 
matrix approach is rather similar to the Born term map, 
due to a limited influence of the non affine component. 

Therefore, our conclusion is that one can safely choose 
among the three approaches to extract correct values for 
the local bulk modulus, which has a small non affine com- 
ponent. However, in the case of the local shear modulus, 
which has a larger non affine component, only the fully 
local approach is appropriate, where both local stress and 
local strain fields are dealt with, and there are no con- 
straints which can limit the correct evolution of the non 
affine component. In this study, we used the fluctuation 
formula Eq. (JSJ to implement the fully local approach. 
An alternative is to obtain the local stress and the lo- 
cal strain relation directly, as was achieved in references 
U [26[ . In this case, one needs to measure, in addition to 
the local stress field crj™, a local strain field e™, derived 
from a coarse-grained displacement field. 
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